﻿
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace LxCFD
{
    public class PressureSIMPLE : ScalarBase
    {
        //local
        public double InitiialPressureProductBoundArea { get; set; }
        public double BoundArea { get; set; }

        public int CenterCellIndex;
        public Vector3d Center;//这个地方如果用get的话，get到的是个副本，会出错
        protected double MinLength;
        //not local
        VelocityVector VelocityVector;
        protected double[] MassFlux;
        Vector3d[] Dc;
        Value1Base<double> Density;
        public void SetZone(Value1Base<double> iniValue, VelocityVector velocityVector, PressureCorrectionBasic Pcor, Value1Base<double> density)
        {
            this.IniValue = iniValue;

            this.VelocityVector = velocityVector;
            this.MassFlux = velocityVector.MassFlux;

            this.Dc = Pcor.Dc;

            this.Density = density;
        }
        public override void Initialize()
        {
            MultiThreads.RunCalTasks(initialize_Setx_CalCenter_step1);
            //计算中几何中心
            {
                ThreadsMath.GetSum(out double sum1, out double sum2, out double sum3);
                Center = new Vector3d(sum1, sum2, sum3);
                Center.MultiplyToSelf(1 / Zone.Volumn);
                MinLength = double.MaxValue;
            }
            MultiThreads.RunCalTasks(initialize_step3);
            //
            {
                ThreadsMath.GetSum(out double sum1, out double sum2, out double sum3);
                InitiialPressureProductBoundArea = sum1;
                BoundArea = 0;
                foreach (var b in Zone.BoundThreads)
                {
                    BoundArea += b.Area;
                }
            }
        }
        protected void initialize_Setx_CalCenter_step1(int taskIndex)
        {
            set_x(IniValue, InnerIndex, taskIndex);
            foreach (var b in VelocityVector.NoSlipWalls)
            {
                set_x(IniValue, b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.Symmetries)
            {
                set_x(IniValue, b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.VelocityInlets)
            {
                set_x(IniValue, b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.PressureInOuts)
            {
                if (b.SetPressureState.GetIsNeedToCal(taskIndex))
                {
                    set_x(b.Pressure, b.ThreadIndex, taskIndex);
                }
            }
            //为了求中心，先求体积*坐标之和
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                Vector3d centerProductVolumn = new Vector3d();
                for (; i1 < endp1; i1++)
                {
                    ref FvCell cell = ref Cells[i1];
                    centerProductVolumn.AddToSelf(cell.Volumn * cell.Center);
                }
                ThreadsMath.AddToSum(centerProductVolumn.X, centerProductVolumn.Y, centerProductVolumn.Z);
            }
        }
        public void CalPressureArea(int taskIndex)     //计算边界压力*面积
        {
            foreach (var b in Zone.BoundThreads)
            {
                b.ThreadIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                double temp = 0;
                for (; i1 < endp1; i1++)
                {
                    ref FvFace face = ref Faces[i1];
                    temp += x[face.NSideCell] * face.Area;
                }
                ThreadsMath.AddToSum(temp);
            }
        }
        public void initialize_step3(int taskIndex)//和上面计算面积的要对应
        {
            CalPressureArea(taskIndex);
            //需要初始化梯度场
            CalGradient(taskIndex);
            //计算找到与几何中心最近的cell
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                double minLengthlocal = double.MaxValue;
                int minIndexLocal = 0;
                for (; i1 < endp1; i1++)
                {
                    double lengthtemp = (Cells[i1].Center - Center).CalLengthSquared();
                    if (lengthtemp < minLengthlocal)
                    {
                        minLengthlocal = lengthtemp;
                        minIndexLocal = i1;
                    }
                }
                lock (MultiThreads.GlobalLock)
                {
                    if (minLengthlocal < MinLength)
                    {
                        MinLength = minLengthlocal;
                        CenterCellIndex = minIndexLocal;
                    }
                }
            }
        }
        public void CorrectPressureByReferenct(double pressureBoundCor, int taskIndex)
        {
            TotalIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                x[i1] += pressureBoundCor;
            }
        }
        public void Predict_bound(int taskIndex)
        {
            foreach (var b in VelocityVector.NoSlipWalls)
            {
                predict_bound_Wall(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.Symmetries)
            {
                predict_bound_Wall(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.VelocityInlets)
            {
                predict_bound_byMassFlux(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.PressureInOuts)
            {
                if (b.SetPressureState.GetIsNeedToCal(taskIndex))
                {
                    set_x(b.Pressure, b.ThreadIndex, taskIndex);
                }
            }
        }
        protected void predict_bound_Wall(IndexNumber indexNumber, int taskIndex)//方法2，精度更高，计算量大，对应P603（有修改，见笔记）
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double fenzi = Vector3d.DotProduct(VelocityVector.GetVelocity(face.OSideCell), face.Normal) * face.OToN.CalLengthSquared();//Dc已改为不包含relax的了
                double fenmu = Vector3d.DotProduct(Vector3d.XXYYZZProduct(Dc[face.OSideCell], face.OToN), face.Normal);
                x[face.NSideCell] = x[face.OSideCell] + Vector3d.DotProduct(Grad[face.OSideCell], face.OToN) + fenzi / fenmu;
            }
        }
        protected void predict_bound_byMassFlux(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double fenzi = (Vector3d.DotProduct(VelocityVector.GetVelocity(face.OSideCell), face.Normal) - MassFlux[i1] / Density.GetFaceValue(i1)) * face.OToN.CalLengthSquared();//Dc已改为不包含relax的了
                double fenmu = Vector3d.DotProduct(Vector3d.XXYYZZProduct(Dc[face.OSideCell], face.OToN), face.Normal);
                x[face.NSideCell] = x[face.OSideCell] + Vector3d.DotProduct(Grad[face.OSideCell], face.OToN) + fenzi / fenmu;
            }
        }



    }
}
//foreach (var b in VelocityVector.NoSlipWalls)
//{
//    BoundArea += Zone.GetThread(b.ThreadID).Area;
//}
//foreach (var b in VelocityVector.Symmetries)
//{
//    BoundArea += Zone.GetThread(b.ThreadID).Area;
//}
//foreach (var b in VelocityVector.VelocityInlets)
//{
//    BoundArea += Zone.GetThread(b.ThreadID).Area;
//}
//foreach (var b in VelocityVector.NoSlipWalls)
//{
//    calPressureArea(b.ThreadIndex, taskIndex);
//}
//foreach (var b in VelocityVector.Symmetries)
//{
//    calPressureArea(b.ThreadIndex, taskIndex);
//}
//foreach (var b in VelocityVector.VelocityInlets)
//{
//    calPressureArea(b.ThreadIndex, taskIndex);
//}
//switch (b.Pressure.valueState)
//{
//    case ValueStates.算一次_未计算:
//        set_x(b.Pressure, b.indexNumber, taskIndex);
//        TasksControl.Sync(taskIndex, b.Pressure.setStateUsed);
//        break;
//    case ValueStates.每次都算:
//        set_x(b.Pressure, b.indexNumber, taskIndex);
//        break;
//    case ValueStates.算一次_已计算:
//    default:
//        break;
//}
//double InitialBoundMeanPressure;
//double CurrentBoundMeadPressure;
//calMeanPressure_Sync(ref InitialBoundMeanPressure, taskIndex);
//double pressureBoundCor = InitialBoundMeanPressure - CurrentBoundMeadPressure;
//protected void calMeanPressure_Sync(ref double pressure, int taskIndex)
//{
//    TasksControl.Sync(taskIndex);
//    foreach (var b in velocityVector.NoSlipWalls)
//    {
//        StatisticTool.AddRange_VolumnWeighted(x, cells, b.indexNumber, taskIndex);
//    }
//    foreach (var b in velocityVector.Symmetries)
//    {
//        StatisticTool.AddRange_VolumnWeighted(x, cells, b.indexNumber, taskIndex);
//    }
//    foreach (var b in velocityVector.VelocityInlets)
//    {
//        StatisticTool.AddRange_VolumnWeighted(x, cells, b.indexNumber, taskIndex);
//    }
//    TasksControl.FinishStatisticAndReset(taskIndex, StatisticTool.MakeResult_WeightedAverage, ref pressure);
//}
//如果边界没有压力参考，取边界平均，然后修正到zone，防止压力跑偏
//if (velocityVector.PressureInOuts.Count == 0)
//{
//    //计算新的平均压力
//    calMeanPressure_Sync(ref CurrentBoundMeadPressure, taskIndex);
//    //并修正
//    double pressureBoundCor = InitialBoundMeanPressure - CurrentBoundMeadPressure;
//    zone.TotalIndex.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        x[i1] += pressureBoundCor;
//    }
//}
//else
//{
//}
//double fenzi = UVW_relaxFactor * Vector3d.DotProduct(velocityVector.getVelocity(face.OwnerCellIndex), face.Normal) * face.OwnerToNeighbor.CalLengthSquared();
//public void calGradPressure(int taskIndex)
//{
//    int i1, endp1;
//    innerIndex.getCellInterval(out i1, out endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        gradPressure[i1] = cells[i1].calCellGradient_LeastSquare(faces, Pressure);
//    }
//    TasksController.syncTasks(taskIndex);
//}